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Abstract 


In  this  paper,  we  present  an  efficient  semi-Lagrangian  based  particle 
level  set  method  for  the  accurate  capturing  of  interfaces.  This  method 
retains  the  robust  topological  properties  of  the  level  set  method  with¬ 
out  the  adverse  effects  of  numerical  dissipation.  Both  the  level  set 
method  and  the  particle  level  set  method  typically  use  high  order  ac¬ 
curate  numerical  discretizations  in  time  and  space,  e.g.  TVD  Runge- 
Kutta  and  HJ-WENO  schemes.  We  demonstrate  that  these  computa¬ 
tionally  expensive  schemes  are  not  required.  Instead,  fast,  low  order 
accurate  numerical  schemes  suffice.  That  is,  the  addition  of  particles  to 
the  level  set  method  not  only  removes  the  difficulties  associated  with 
numerical  diffusion,  but  also  alleviates  the  need  for  computationally 
expensive  high  order  accurate  schemes.  We  use  an  efficient,  first  order 
accurate  semi-Lagrangian  advection  scheme  coupled  with  a  first  order 
accurate  fast  marching  method  to  evolve  the  level  set  function.  To 
accurately  track  the  underlying  flow  characteristics,  the  particles  are 
evolved  with  a  second  order  accurate  method.  Since  we  avoid  complex 
high  order  accurate  numerical  methods,  extending  the  algorithm  to 
arbitrary  data  structures  becomes  more  feasible,  and  we  show  prelim¬ 
inary  results  obtained  with  an  octree-based  adaptive  mesh. 

‘Research  supported  in  part  by  an  ONR  YIP  and  PECASE  award  (N00014-01-1-0620), 
a  Packard  Foundation  Fellowship,  a  Sloan  Research  Fellowship,  ONR  N00014-03- 1-0071, 
ONR  N00014-02-1-0720,  NSF  DMS-0106694,  NSF  ITR-0121288  and  DOE  under  the  ASCI 
Academic  Strategic  Alliances  Program  (LLNL  contract  B341491).  In  addition,  the  first 
author  was  supported  in  part  by  an  NSF  postdoctoral  fellowship  (DMS-0202459). 

^Mathematics  Department,  UCLA,  Los  Angeles,  CA  90095. 

^Computer  Science  Department,  Stanford  University,  Stanford,  CA  94305. 


1 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

25  APR  2004 


2.  REPORT  TYPE 


4.  TITLE  AND  SUBTITLE 

A  Fast  and  Accurate  Semi-Lagrangian  Particle  Level  Set  Method 


6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Stanford  University  , Computer  Science  Department, Stanford, CA, 94305 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED 

00-00-2004  to  00-00-2004 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

In  this  paper,  we  present  an  efficient  semi-Lagrangian  based  particle  level  set  method  for  the  accurate 
capturing  of  interfaces.  This  method  retains  the  robust  topological  properties  of  the  level  set  method 
without  the  adverse  effects  of  numerical  dissipation.  Both  the  level  set  method  and  the  particle  level  set 
method  typically  use  high  order  ac-  curate  numerical  discretizations  in  time  and  space,  e.g.  TVD 
Runge-Kutta  and  HJ-WENO  schemes.  We  demonstrate  that  these  computationally  expensive  schemes  are 
not  required.  Instead,  fast,  low  order  accurate  numerical  schemes  suffice.  That  is,  the  addition  of  particles 
to  the  level  set  method  not  only  removes  the  difficulties  associated  with  numerical  diffusion,  but  also 
alleviates  the  need  for  computationally  expensive  high  order  accurate  schemes.  We  use  an  efficient,  first 
order  accurate  semi-Lagrangian  advection  scheme  coupled  with  a  first  order  accurate  fast  marching 
method  to  evolve  the  level  set  function.  To  accurately  track  the  underlying  flow  characteristics,  the 
particles  are  evolved  with  a  second  order  accurate  method.  Since  we  avoid  complex  high  order  accurate 
numerical  methods,  extending  the  algorithm  to  arbitrary  data  structures  becomes  more  feasible,  and  we 
show  preliminary  results  obtained  with  an  octree-based  adaptive  mesh. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

24 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


1  Introduction 


Standard  Eulerian  advection  algorithms,  e.g.  HJ-(W)ENO  methods  [18,  8] 
combined  with  TVD-based  high  order  accurate  Runge-Kutta  schemes  [18], 
require  a  strict  bound  on  the  maximum  possible  time  step  due  to  a  stability- 
based  CFL  criterion.  On  the  other  hand,  particles  are  not  restricted  by  this 
stability  criterion,  rather  the  size  of  the  time  step  used  can  be  based  solely 
on  the  degree  of  numerical  accuracy  desired.  Grid  based  semi-Lagrangian 
advection  methods  are  likewise  not  limited  by  a  stability-based  CFL  condi¬ 
tion,  since  each  grid  point  is  treated  in  a  “particle-like”  manner.  However, 
semi-Lagrangian  schemes  can  suffer  from  large  amounts  of  numerical  dis¬ 
sipation,  making  their  use  problematic.  HJ-(W)ENO  methods  experience 
far  less  numerical  dissipation  due  to  their  high  order  accurate  adaptive  na¬ 
ture.  In  order  to  take  advantage  of  the  stability  afforded  by  particle  based 
methods,  the  spatial  and  temporal  coherency  of  Eulerian  methods,  and  the 
opportunity  for  selective  adaptive  mesh  refinement  near  the  interface  in  or¬ 
der  to  resolve  small  scale  features  as  discussed  below,  we  propose  to  couple 
together  a  semi-Lagrangian  advection  scheme  with  a  characteristic-based 
particle  method  to  track  a  passively  advected  interface. 

A  flexible  and  easy-to-implement  interface  tracking  technique  is  the  level 
set  method  of  Osher  and  Sethian  [12].  By  storing  the  distance  to  the  in¬ 
terface  at  each  point  on  a  fixed  computational  grid,  handling  gross  changes 
to  interface  topology,  e.g.  pinching  and  merging,  becomes  trivial  as  com¬ 
pared  to  standard  Lagrangian  techniques  [25]  that  typically  require  ad  hoc 
techniques  to  address  mesh  connectivity  during  merging  and  pinching.  By 
avoiding  these  difficulties  and  utilizing  well  established  numerical  algorithms 
for  the  solution  of  nonlinear  hyperbolic  conservation  laws,  level  set  methods 
have  been  applied  to  a  wide  variety  of  problems  including  fluid  mechanics, 
computer  vision,  material  science,  and  computer  graphics.  One  difficulty 
with  the  use  of  the  level  set  method  is  the  need  to  control  the  numerical 
diffusion  (or  mass  loss)  present  in  the  method,  especially  in  areas  of  high 
curvature  and  long,  thin  filamentary  regions.  Various  authors  [24,  22,  23] 
have  attempted  to  correct  this  problem  by  reinitializing  the  level  set  func¬ 
tion,  4>,  to  be  signed  distance  to  the  interface  after  each  time  step.  High 
order  accurate  TVD  Runge-Kutta  and  HJ-(W)ENO  techniques  can  be  used 
to  perform  this  reinitialization.  While  producing  reasonable  results,  these 
methods  suffer  from  the  same  Eulerian-based  advection  issues  mentioned 
above  such  as  small  time  step  restrictions. 

An  alternative  characteristic-based  diffusion  correction  technique,  the 
particle  level  set  method  [5],  has  been  recently  proposed.  In  this  method, 
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two  sets  of  marker  particles  are  placed  near  the  interface,  one  set  associated 
with  the  interior  (</>  <  0)  region,  and  the  other  with  the  exterior  (c j>  >  0) 
region.  Errors  due  to  numerical  dissipation  can  then  be  identified  when 
interior  particles  appear  in  the  exterior  region  or  exterior  particles  appear 
in  the  interior  region.  Since  the  particles  are  able  to  more  accurately  track 
the  underlying  flow  characteristics,  these  “escaped”  particles  can  be  used  to 
correct  the  level  set  representation  of  the  interface.  The  particle  level  set 
method  has  been  shown  to  possess  excellent  volume  conservation  properties 
and  a  high  degree  of  geometrical  accuracy  in  tracking  contact  discontinuities, 
comparable  to  other  interface  methods  including  Volume-Of-Fluid  (VOF) 
and  explicit  front  tracking.  At  the  same  time  the  method  maintains  the 
highly  desirable  topological  properties  and  ease-of-implementation  of  the 
original  level  set  method.  As  an  example  of  the  flexibility  of  the  particle  level 
set  method,  its  use  in  modeling  complex  three  dimensional  water  surfaces 
can  be  seen  in  [6,  7]. 

In  [5],  HJ-WENO  numerical  methods  were  used  to  evolve  and  reinitialize 
(f>.  Reinitialization  of  4>  was  used  in  order  to  assist  the  particles  in  obtain¬ 
ing  an  accurate  distance  to  the  interface.  Due  to  the  HJ-WENO  advection 
scheme  used,  a  stability  based  CFL  condition  was  imposed  on  the  size  of 
the  time  step.  Use  of  this  scheme  in  an  adaptive  mesh  setting  would  im¬ 
pose  a  severe  time  step  restriction  on  a  simulation,  discouraging  the  use  of 
selective  mesh  adaptation  to  resolve  small  scale  features.  Due  to  the  mini¬ 
mal  amount  of  diffusion  exhibited  by  the  particle  level  set  method,  we  can 
avoid  the  common  but  computationally  costly  approach  of  adapatively  re¬ 
solving  the  mesh  just  for  the  sake  of  minimizing  numerical  diffusion  at  the 
interface.  Moreover,  since  the  particles  dictate  a  sharp  and  geometrically 
accurate  interface,  we  propose  to  use  a  first  order  accurate  semi-Lagrangian 
advection  method  [4,  15,  19].  Despite  being  unconditionally  stable,  use  of  a 
first  order  accurate  sem-Lagrangian  scheme  has  typically  been  shunned  due 
to  the  large  amount  of  numerical  diffusion  inherently  incurred.  Although  a 
low  order  accurate  semi-Lagrangian  scheme  was  used  by  [20]  for  level  set  ad¬ 
vection,  computationally  expensive  higher  order  accurate  semi-Lagrangian 
methods  are  usually  preferred,  see  [21].  Due  to  the  observed  excellent  dif¬ 
fusion  limiting  properties  of  the  particles,  we  illustrate  that  the  fast  first 
order  accurate  semi-Lagrangian  scheme  is  sufficient.  Also,  we  replace  the 
HJ-WENO  reinitialization  scheme  with  a  fast,  i.e.  0(N  log  IV),  first  order 
accurate  marching  technique  first  proposed  by  Tsitsiklis  [26]  and  later  pop¬ 
ularized  by  Sethian  and  co-workers,  see  e.g.  [17]. 

None  of  the  methods  we  propose  are  bound  by  a  grid  based  CFL  stability 
condition,  rather  only  numerical  accuracy  needs  to  be  taken  into  account. 
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The  resulting  numerical  method  is  a  computationally  fast  and  geometrically 
accurate  interface  tracking  technique  that  efficiently  provides  for  the  adap¬ 
tive  resolution  of  small  scale  features.  We  demonstrate  this  last  claim  by 
showing  some  preliminary  computations  on  an  octree  data  structure. 

2  Numerical  Method 

2.1  Level  Set  Method 

The  underlying  idea  behind  level  set  methods  is  to  embed  an  interface  T 
which  bounds  a  region  O  C  R3  as  the  zero  level  set  of  a  higher  dimensional 
function  cf>(x,t).  The  level  set  function  has  the  following  properties, 

4>(x,  t)  >  0  for  x  0  n 
(f)(x,  t)  <  0  for  x  €  12, 

where  we  include  (f>  =  0  with  the  negative  4>  values.  Then  the  interface  lies 
in  between  cj>  >  0  and  <f>  =  0,  but  can  of  course  be  identified  as  0  =  0. 
Note  that  (/>  is  a  scalar  function  in  R3  which  greatly  reduces  the  complexity 
of  describing  the  interface,  especially  when  undergoing  topological  changes 
such  as  pinching  and  merging. 

The  motion  of  the  interface  is  determined  by  a  velocity  field,  u,  which 
can  depend  on  a  variety  of  things  including  position,  time,  geometry  of  the 
interface,  or  be  given  externally  for  instance  as  the  material  velocity  in  a 
fluid  flow  simulation.  In  most  of  the  examples  below,  the  velocity  field  is 
externally  given,  and  the  evolution  equation  for  the  level  set  function  is  given 

by 

4>t  +  u  ■  V(f>  =  0.  (1) 

In  order  to  allow  for  a  computationally  efficient  implementation,  we  solve 
this  equation  locally  near  the  interface  in  a  manner  similar  to  [1,  13].  We 
solve  equation  (1)  in  a  region  of  ±5  max(Ax,  Ay)  of  the  interface.  The  size 
of  this  region  is  chosen  such  that  the  (f>  =  0  level  set  is  not  advected  outside 
this  region  in  one  semi-Lagrangian  time  step  as  discussed  below. 

It  is  convenient  to  initialize  ^  to  be  a  signed  distance  function  with 
| V0|  =  1.  This  ensures  that  the  level  set  is  a  smoothly  varying  function 
well  suited  for  accurate  numerical  computations.  Unfortunately,  as  noted  in 
[24],  the  level  set  function  can  quickly  cease  to  be  a  signed  distance  function 
especially  for  flows  undergoing  extreme  topological  changes.  Reinitialization 
algorithms  maintain  the  signed  distance  property  by  solving  to  steady  state 
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(as  fictitious  time  r  — >  oo)  the  equation 


4>t  +  sgn(^o)(|V0|  -  1)  =  0  (2) 


where  sgn(c/>o)  is  a  one-dimensional  smeared  out  signurn  function  approxi¬ 
mated  numerically  in  [24]  as 


sgn(0o) 


<t>o 

V^o  +  (Ax’)2 


Efficient  ways  to  solve  equation  (2)  to  steady  state  via  fast  marching  methods 
are  discussed  in  [16,  17].  Again,  equation  (2)  only  needs  to  be  solved  locally 
near  the  interface.  We  reinitialize  the  level  set  function  via  equation  (2)  in 
the  same  region  about  the  interface  as  we  solve  equation  (1)  in. 

Geometric  quantities  are  easily  calculated  using  the  level  set  function, 
including  the  unit  normal, 


(3) 


and  the  curvature, 


k  =  V  • 


(4) 


The  spatial  derivatives  in  equations  (3)  and  (4)  can  be  calculated  using 
standard  central  differencing  operators  when  the  denominators  are  non-zero. 
Otherwise,  one  sided  differencing  is  used.  For  more  details  on  level  set 
methods,  we  refer  the  interested  reader  to  [11]. 


2.2  Particle  Level  Set  Method 

The  particle  level  set  method  [5]  is  a  thickened  front  tracking  approach 
which  uses  particles  to  assist  the  level  set  method  in  accurately  tracking  flow 
characteristics  in  under-resolved  regions  (and  consequently  preserve  mass). 
This  is  achieved  through  the  placement  of  massless  marker  particles  near 
the  interface  as  a  diffusion  correction  mechanism  for  the  level  set  function. 

Two  sets  of  marker  particles  are  randomly  placed  in  a  “thickened”  sur¬ 
face  region  about  the  (f)  =  0  level  set.  The  thickness  of  the  band  used  in  the 
examples  section  is  three  grid  cells  on  each  side  of  the  interface.  Positive 
particles  are  located  in  the  (f>  >  0  region  and  negative  particles  in  the  <j>  <  0 
region.  The  number  of  particles  placed  in  each  cell  can  be  adjusted  accord¬ 
ing  to  the  amount  of  surface  resolution  desired.  For  the  examples  presented, 
16  particles  per  cell  were  used  in  2D  and  32  in  3D  as  suggested  in  [5] .  Each 
particle  possesses  a  radius,  rp.  which  is  constrained  between  a  minimum  and 
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maximum  value  based  upon  the  size  of  the  underlying  computational  grid.  A 
minimum  radius  of  .lmin(Ax,  Ay)  and  maximum  radius  of  .5  min(  Ax,  Ay) 
is  used.  The  radius  of  a  particle  changes  dynamically  throughout  the  sim¬ 
ulation  since  the  particle’s  relative  location  to  the  surface  changes  in  time. 
This  allows  for  a  multiscale  sampling  of  the  interface  by  the  particles.  The 
radius  of  each  particle  is  set  according  to: 

max  if  Sp(/)(xp)  > 

V  max 

Sp<j>(Xp)  if  V min  ^  ^  ^ max  (5) 

^  min  if  sp(j)(xp)  < 

T  mini 

where  sp  is  the  sign  of  the  particle  (+1  for  positive  particles  and  -1  for  neg¬ 
ative  particles).  This  radius  adjustment  keeps  the  boundary  of  the  particles 
tangent  to  the  surface  whenever  possible. 


2.2.1  Semi-Lagrangian  Advection 

In  the  original  particle  level  set  method,  high  order  accurate  numerical 
schemes,  i.e.  3rd  order  accurate  TVD  Runge-Kutta  [18]  in  time  and  5th 
order  accurate  Hamilton- Jacobi  WENO  [8]  for  the  advection  term,  were 
used  to  evolve  the  level  set  function.  Also,  the  particles  were  integrated 
forward  in  time  with  a  3rd  order  accurate  TVD  Runge-Kutta  method  with 
bilinear  interpolation  used  to  calculate  the  velocity  of  a  given  particle  on  the 
computational  grid.  We  replace  the  high  order  advection  and  time  integra¬ 
tion  schemes  with  a  fast  first  order  accurate  semi-Lagrangian  method,  e.g. 
given  a  computational  grid,  Xij  =  ( iAx,jAy ),  and  temporal  discretization, 
tn  =  nAt, 


(fiij1  —  a/3</>”+ls+1  +  (l  — a)/?</>”s+1  +  a(l  — /3)(/>”+ls  +  (l  — a)(l  — (6) 
where 


r  =  i  — 


s  =  3 


At 

Ui’jA x 

At 
Ay 


■UJ 


a  = 


_  (i-r)Ax-UijAt 


Ax  ’ 


P  = 


Ay 


(7) 

(8) 


and  u(xi,j )  =  (uij,Vij).  This  method  is  unconditionally  stable  due  to  the 
linear  interpolation  of  <f>  in  equation  (6)  and  is  convergent  to  the  correct 
solution  according  to  the  Lax-Richtmyer  theorem.  Further  discussion  of 
this  scheme  can  be  found  in  [20].  Due  to  the  unconditional  stability  of  this 
scheme,  the  size  of  the  time  step  is  not  governed  by  a  stability-based  CFL 
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condition.  A  CFL  number  of  4.9  (i.e.  Af(|u|maa;/Ax  +  \v\max/A y)  =  4.9) 
was  used  in  the  examples  section. 

Adequate  resolution  of  the  underlying  flow  field  by  the  particles  is  nec¬ 
essary  in  order  to  maintain  an  accurate  representation  of  the  interface.  We 
have  found  that  a  second  order  accurate  Runge-Kutta  midpoint  rule  is  re¬ 
quired  for  the  time  integration  of  the  particles.  In  fact,  lowering  the  time 
integration  of  the  particles  from  second  order  to  first  order  does  have  an  ad¬ 
verse  effect  on  the  numerical  results,  even  though  replacing  the  high  order 
accurate  integration  of  the  level  set  equation  with  the  first  order  accurate 
semi-Lagrangian  and  fast  marching  methods  seems  to  have  little  effect.  Bi¬ 
linear  interpolation  is  used  to  calculate  the  velocity  of  a  particle  from  the 
computational  grid.  The  time  step  used  for  the  particles  is  the  same  as  that 
used  for  the  level  set  advection. 

2.2.2  Error  Correction 

After  the  level  set  function,  and  both  the  positive  and  negative  particles 
are  integrated  separately  forward  in  time,  the  particles  are  used  to  correct 
any  errors  in  the  representation  of  the  interface  according  to  the  level  set 
function.  This  particle  correction  mechanism  is  comprised  of  several  steps 
discussed  below.  Note  that  we  apply  the  error  correction  step  after  each 
modification  of  the  level  set,  i.e.  after  the  advection  step  and  after  the  reini¬ 
tialization  step. 

Identification  of  Error  When  particles  appear  on  the  wrong  side  of 
the  interface  by  more  than  their  radius,  we  indicate  the  presence  of  an  error 
in  the  level  set  representation  of  the  interface.  These  particles  are  said  to 
have  escaped.  In  smooth,  well  resolved  regions  of  the  flow  where  the  level  set 
method  is  accurate,  the  particles  do  not  drift  an  appreciable  amount  across 
the  interface,  so  we  do  not  use  the  particle  representation  of  the  interface  in 
this  region.  Only  when  the  semi-Lagrangian  advection  of  the  level  set  has 
clearly  made  an  error  do  we  resort  to  using  the  following  steps  to  reconstruct 
the  level  set  function. 

Quantification  of  Error  For  each  particle  p.  we  associate  a  spherical 
level  set  function,  (f>p ,  whose  size  is  determined  by  the  particle  radius,  i.e. 

(t>p{x)  =  sp(rp  -  \x  -  xp\).  (9) 

The  particle  defined  level  set  function  is  computed  locally  on  the  eight  cor¬ 
ners  of  the  cell  containing  the  particle.  The  local  values  of  (f)p  are  the  particle 
predictions  of  the  values  of  the  overall  level  set  function,  <f>,  on  the  corners 
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of  the  cell.  Any  variation  of  <fi  from  c j>p  indicates  potential  errors  in  the  level 
set  solution. 

Error  Correction  We  use  the  escaped  positive  particles  to  rebuild  the 
(f>  >  0  region  and  the  escaped  negative  particles  to  rebuild  the  (f>  <  0  region. 
For  example,  take  the  (f>  >  0  region  and  an  escaped  positive  particle.  Using 
equation  (9),  the  (j)p  values  of  the  eight  grid  points  on  the  boundary  of  the 
cell  containing  the  particle  are  calculated.  Each  4>p  is  compared  to  the  local 
value  of  (j)  and  the  maximum  of  these  two  values  is  taken  as  4>+ .  This  is  done 
for  all  escaped  positive  particles,  creating  a  reduced  error  representation  of 
the  4>  >  0  region.  That  is,  given  a  level  set  4>  and  a  set  of  escaped  positive 
particles  E+,  we  initialize  (f>+  with  cj)  and  then  calculate 

4>+ =  max  (f>+).  (10) 

Vp€E+ 

Similarly,  to  calculate  a  reduced  error  representation  of  the  (f>  <  0  region, 
we  initialize  <\f~  with  <fi  and  then  calculate 


(j)  =  min  ). 


(11) 


< f>+  and  4>-  will  not  agree  due  to  the  errors  in  both  the  particle  and  level  set 
methods  as  well  as  interpolation  errors,  etc.  We  merge  (f>+  and  (f>~  back  into 
a  single  level  set  by  setting  cj)  equal  to  the  value  of  (f>+  or  (j)~  which  is  least 
in  magnitude  at  each  grid  point, 


/  </>+  if  |</>+|  <  I <t>~ 
\  4>~  if  l</>+l  >  1 4>~ 


(12) 


The  minimum  magnitude  is  used  to  reconstruct  the  interface  (instead  of,  for 
example,  taking  an  average),  since  it  gives  priority  to  values  that  are  closer 
to  the  interface. 


2.2.3  Reinitialization  and  Radii  Adjustment 

(j)  is  maintained  to  be  a  signed  distance  function  by  solving  equation  (2)  via 
a  fast  marching  technique.  Again,  for  the  sake  of  efficiency,  we  only  reini¬ 
tialize  cj)  within  a  band  of  the  interface.  Combining  this  narrow  banding 
optimization  with  the  fast  marching  method  provides  for  a  very  fast  reini¬ 
tialization  procedure.  To  ensure  proper  <j>  values  for  the  semi-Lagrangian 
update,  we  reinitialize  (j)  within  a  band  of  ±6max(Ax,  A y)  of  the  interface. 
This  procedure  is  performed  after  each  combined  semi-Lagrangian  update 
and  error  correction  step.  Unfortunately,  reinitialization  may  cause  the  zero 


level  set  to  move,  which  is  not  desirable,  so  we  use  the  particles  to  correct 
these  errors  as  well.  Finally,  the  particles  resample  their  position  relative  to 
the  0  =  0  level  set  and  adjust  their  radii  accordingly.  Any  particles  which 
remain  escaped  have  their  radius  set  to  the  minimum  particle  radius  value. 

In  summary,  the  order  of  operations  is:  evolve  both  the  particles  and  the 
level  set  function  forward  in  time,  correct  errors  in  the  level  set  function  using 
particles,  apply  the  fast  marching  method  to  reinitializae  0  in  a  band  near 
the  interface,  again  correct  errors  in  the  level  set  function  using  particles, 
and  finally  adjust  the  particle  radii. 

3  Examples 

3.1  Rigid  Body  Rotation  of  Zalesak’s  Disk 

Consider  the  rigid  body  rotation  of  Zalesak’s  disk  in  a  constant  vorticity 
velocity  field  [27].  The  initial  data  is  a  slotted  circle  centered  at  (50,75) 
with  a  radius  of  15,  a  width  of  5,  and  a  slot  length  of  25.  The  constant 
vorticity  velocity  field  is  given  by 

u  =  (vr/314)  (50  —  y), 

v  =  (7r/314)(x  —  50), 

so  that  the  disk  completes  one  revolution  every  628  time  units. 

To  better  understand  the  ability  of  the  various  coupled  advection  and 
reinitialization  algorithms  to  accurately  respresent  a  passively  advected  in¬ 
terface,  a  comparison  of  a  level  set  only  method  utilizing  these  algorithms 
has  been  performed.  The  use  of  the  expensive,  but  accurate  HJ-WENO 
schemes  in  a  level  set  only  method  is  clearly  justified  by  the  results  seen  in 
table  1  and  figure  1.  However,  the  ability  of  an  HJ-WENO  advection  scheme 
to  accurately  advect  the  interface  can  be  severely  impaired  when  coupled 
with  a  low  order  accurate  reinitialization  method  such  as  the  fast  marching 
method  as  seen  in  figure  1(b).  The  low  order  accurate  errors  introduced  by 
the  use  of  a  fast  marching  reinitialization  method  are  compounded  due  to  the 
small  time  step  used  during  the  HJ-WENO  advection  phase,  reducing  the 
overall  accuracy  of  the  coupled  method.  Semi-Lagrangian  based  methods 
minimize  this  effect  due  to  the  larger  time  step  allowed.  That  is,  the  larger 
time  step  means  that  the  diffusion  errors  from  the  spatial  mis-approximation 
of  the  interface  are  applied  less  often. 

Figures  2  and  3  compare  the  evolution  of  a  high  order  accurate  level 
set  only  method  (3rd  order  TYD  RK  in  time  and  5th  order  HJ-WENO  in 
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Advection  Method  + 
Reinitialization  Method 

Area 

%  Area  Loss 

CPU  Time  (sec) 

exact 

582.2 

- 

- 

HJ-WENO  +  HJ-WENO 

604.4 

-4.08% 

86.36 

HJ-WENO  +  FMM 

181.3 

68.78% 

11.39 

SL  +  HJ-WENO 

328.7 

43.40% 

7.20 

SL  +  FMM 

298.8 

48.55% 

.45 

Table  1:  Zalesak’s  disk.  Level  set  only  method  comparing  different  advection 
and  reinitialization  schemes  on  a  100x100  computational  cell  grid  after  one 
rotation. 

space),  the  original  high  order  accurate  particle  level  set  method,  and  the 
newly  proposed  fast  semi-Lagrangian  particle  level  set  method  (coupled  to 
a  fast  marching  method),  after  one  and  two  revolutions,  respectively.  The 
exact  solution  is  also  plotted  for  the  sake  of  comparison.  As  expected,  the 
level  set  only  method  applies  an  excessive  amount  of  regularization  in  the 
sharp  corners  which  the  particles  correct.  The  ability  of  the  particles  to 
maintain  sharp  features,  even  when  a  highly  diffusive  first  order  accurate 
semi-Lagrangian  advection  scheme  is  combined  with  a  first-order  accurate 
fast  marching  reinitialization  method,  is  quite  remarkable.  Again,  we  note 
that  a  second  order  Runge  Kutta  midpoint  scheme  is  used  for  the  particle 
advection.  Use  of  a  higher-order  accurate  time  integration  scheme  for  the 
particles  or  a  HJ-WENO  based  advection  and/or  reintialization  scheme  did 
not  significantly  add  to  the  quality  of  the  solution  already  obtained. 

Tables  2  and  3  compare  the  area  loss  (or  gain)  of  the  original  high-order 
accurate  particle  level  set  method  to  the  newly  proposed  semi-Lagrangian 
based  particle  level  set  method  on  three  different  grids.  A  CFL  number  of  4.9 
is  used  on  all  grids  for  the  semi-Lagrangian  calculation,  while  a  CFL  number 
of  .5  is  used  in  all  HJ-WENO  based  calculations.  The  area  is  calculated 
using  a  second  order  accurate  unbiased  level  set  contouring  algorithm  [3]. 
In  addition,  we  calculate  the  accuracy  of  the  interface  location  using  the 
first  order  accurate  error  measure  introduced  in  [22], 

J  I H \(f) expected)  ^  computed)  \  dxdy ,  (13) 

where  L  is  the  length  of  the  expected  interface.  This  integral  is  numerically 
calculated  as  in  [22]: 

•  partition  the  domain  into  many  tiny  pieces  (1000  x  1000), 
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Grid  Cells 

Area 

%  Area  Loss 

L\  Error 

Order 

CPU  Time  (sec) 

exact 

582.2 

- 

- 

- 

- 

One 

Revolution 

50 

562.2 

3.6  % 

.302 

N/A 

25.076 

100 

578.4 

.41% 

.073 

2.04 

134.043 

200 

581.7 

.08% 

.031 

1.22 

840.038 

Two 

Revolutions 

50 

552.2 

5.32% 

.357 

N/A 

50.052 

100 

576.3 

.78% 

.092 

1.95 

267.214 

200 

580.6 

.20% 

.037 

1.33 

1691.55 

Table  2:  Zalesak’s  disk.  HJ-WENO  based  Particle  Level  Set  method. 


Grid  Cells 

Area 

%  Area  Loss 

L\  Error 

Order 

CPU  Time  (sec) 

exact 

582.2 

- 

- 

- 

- 

One 

Revolution 

50 

565.2 

3.09% 

.434 

N/A 

.891 

100 

574.5 

1.07% 

.181 

1.25 

4.045 

200 

580.5 

.22% 

.105 

.79 

18.215 

Two 

Revolutions 

50 

567.7 

2.66% 

.610 

N/A 

1.752 

100 

574.9 

1.01% 

.206 

1.57 

7.941 

200 

580.5 

.22% 

.103 

1.00 

36.341 

Table  3:  Zalesak’s  disk.  Semi-Lagrangian  based  Particle  Level  Set  method. 


•  interpolate  < ^computed  onto  the  newly  partitioned  domain  and  calculate 
4> expected  for  the  domain, 

•  numerically  integrate  equation  (13),  where  H(<f))  is  the  indicator  func¬ 
tion  for  (j)  <  0,  i.e.  H(cj))  =  1  if  </>  <  0  and  H(cj))  =  0  otherwise. 

We  note  that  both  schemes  are  comparable  in  the  quality  of  interface  re¬ 
construction,  while  the  fast  semi-Lagrangian  based  method  is  far  superior 
in  CPU  time  used. 

3.2  Single  Vortex 

While  Zalesak’s  disk  is  a  good  indicator  of  diffusion  errors  in  an  interface 
capturing  method,  it  does  not  test  the  ability  of  an  Eulerian  scheme  to 
accurately  resolve  thin  filaments  on  the  scale  of  the  mesh  which  can  occur 
in  stretching  and  tearing  flows.  A  flow  which  exhibits  interface  stretching  is 
the  “vortex- in-a-box”  problem  introduced  in  [2] .  The  velocity  field  is  defined 
by  the  stream  function 

'L  =  —  sin2(7nc)  sin2(7ry). 

7 r 


11 


Grid  Cells 

Area 

%  Area  Loss 

L\  Error 

Order 

CPU  Time  (sec) 

exact 

.0707 

- 

- 

- 

- 

64 

.0694 

1.68% 

2.89e-3 

N/A 

152.448 

128 

.0701 

.79% 

1.40e-3 

1.0 

885.974 

256 

.0705 

.32% 

5.43e-4 

1.4 

6609.78 

Table  4:  One  period  of  vortex  flow.  HJ-WENO  Particle  Level  Set  method. 


Grid  Cells 

Area 

%  Area  Loss 

L\  Error 

Order 

CPU  Time  (sec) 

exact 

.0707 

- 

- 

- 

- 

64 

.0693 

1.83% 

3.34e-3 

N/A 

4.316 

128 

.0701 

.73% 

9.73e-4 

1.8 

23.653 

256 

.0704 

.38% 

5.58e-4 

.8 

149.995 

Table  5:  One  period  of  vortex  flow.  Semi-Lagrangian  Particle  Level  Set 
method. 


A  unit  computational  domain  is  used  with  a  circle  of  radius  .15  placed  at 
(.5,  .75).  The  resulting  velocity  field  stretches  out  the  circle  into  a  long,  thin 
filament  which  progressively  wraps  itself  towards  the  center  of  the  box.  In 
under-resolved  regions,  the  particles  will  not  be  close  enough  together  to 
accurately  represent  the  interface  and  thin  filament  structures  will  break 
apart.  However,  the  particles  still  track  the  interface  motion  with  second 
order  accuracy,  and  thus  the  resulting  pieces  are  in  accurate  locations. 

For  the  purposes  of  error  analysis,  the  velocity  field  is  time  reversed  by 
multiplying  by  cos(7r t/T)  where  T  is  the  time  at  which  the  flow  returns  to 
its  initial  state,  see  [9].  The  reversal  period  used  in  the  error  analysis  of  the 
vortex  problem  is  T  =  8  producing  a  maximal  stretching  as  seen  in  figure  4. 
As  can  be  seen  from  the  error  tables  4  and  5  as  well  as  figures  5  and  6, 
the  ability  of  the  fast,  first  order  accurate  semi-Lagrangian  particle  level  set 
method  to  model  interfaces  undergoing  substantial  stretching  is  comparable 
to  the  significantly  slower  HJ-WENO  method.  The  L\  errors  reported  here 
for  both  cases  compare  favorably  with  those  reported  in  [14]  using  a  VOF 
PLIC  method.  Again,  the  semi-Lagrangian  method  is  substantially  faster 
due  to  the  much  larger  characteristic-based  CFL  number,  4.9,  used  as  com¬ 
pared  to  the  “safe”  stability-based  CFL  number  of  .5  for  the  HJ-WENO 
method. 
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3.3  Octree  Example 

[9]  proposed  a  three  dimensional  incompressible  flow  field  which  combines  a 
deformation  in  the  x-y  plane  with  one  in  the  x-z  plane.  The  velocity  field  is 
given  by 


u(x,y,z)  =  2  sin2(7rx)  sin(27ry)  sin(27rz), 

v(x,y,z )  =  —  sin(27r.T)  sin2(7ry)  sin(27rz), 

w(x,y,z)  =  —  sin(27r.T)  sin(27ry)  sin2(7T2;) 

and  the  flow  field  is  modulated  in  time  with  a  period  of  T  =  3.  A  sphere  of 
radius  .15  is  placed  within  a  unit  computational  domain  at  (.35,  .35,  .35).  [5] 
used  this  3D  deformation  field  test  to  demonstrate  the  ability  of  the  particles 
to  help  conserve  the  volume  of  the  sphere  as  shown  in  figure  7.  A  uniform 
1003  grid  cell  domain  was  used  in  this  test  case.  However,  this  test  case  also 
shows  some  surface  aliasing  when  the  resolution  of  the  computational  grid 
is  insufficient  to  resolve  the  small  scale  features  and  thin  filaments  formed 
during  the  simulation.  Increasing  the  resolution  of  the  uniform  grid  used  to 
store  cf)  in  order  to  accurately  resolve  the  surface  is  prohibitively  expensive 
in  terms  of  computational  time  and  memory  usage.  Since  particles  track  the 
surface  well  even  when  the  level  set  is  underresolved,  the  original  particle 
level  set  method  is  able  to  show  excellent  volume  conservation  properties  on 
the  coarse  grid  used. 

Octrees  allow  for  significantly  higher  effective  grid  resolutions  and  selec¬ 
tive  parts  of  the  computational  domain  can  be  refined  to  accurately  resolve 
the  surface  where  needed.  Figure  8  shows  the  same  test  using  an  octree 
grid  refinement  scheme.  A  maximum  refinement  depth  of  10,  resulting  in 
an  effective  resolution  of  5123  grid  cells,  was  used.  Use  of  octree  grid  refine¬ 
ment  guarantees  that  the  level  set  is  never  underresolved  anytime  during  the 
simulation  and  that  the  aliasing  problems  that  occur  in  the  lower  resolution 
uniform  grid  simulation  are  non-existent.  The  time  and  memory  required 
for  the  level  set  in  this  example  is  significantly  less  than  running  the  same 
simulation  with  a  highly  refined  uniform  grid. 

We  stress  that  the  octree  based  implementation  and  efficiency  is  facili¬ 
tated  since  we  only  require  a  semi-Lagrangian  advection  scheme  and  a  fast 
marching  method,  i.e.  as  opposed  to  TVD  RK  for  time  evolution  and  HJ- 
WENO  for  both  the  advection  and  the  reinitialization  equation.  For  a  more 
complete  description  of  the  octree  particle  level  set  method,  see  [10]. 


13 


(a)  HJ-WENO  advection  + 
HJ-WENO  reinitializaiton 


(b)  HJ-WENO  advection  + 
FMM  reinitialization 


(c)  Semi-Lagrangian  advection 
+  HJ-WENO  reinitialization 


(d)  Semi-Lagrangian  advection 
+  FMM  reinitialization 


Figure  1:  Level  set  only  method.  Comparison  of  various  advection  and 
reinitialization  algorithms  after  one  revolution  on  a  100  x  100  grid  cell  com¬ 
putational  domain. 
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(c)  HJ-WENO  particle  level  set  (d)  Semi-Lagrangian  particle 

level  set 


Figure  2:  Comparison  of  particle  level  set  methods  after  one  revolution  on 
a  100  x  100  grid  cell  computational  domain. 
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(a)  Initial  notched  disk 


(b)  Level  set 


(c)  HJ-WENO  particle  level  set  (d)  Semi-Lagrangian  particle 

level  set 

Figure  3:  Comparison  after  two  revolutions  on  a  100  x  100  grid  cell  compu¬ 
tational  domain. 
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(a)  HJ-WENO  particle  level  set  (b)  Semi-Lagrangian  particle 

level  set 

Figure  4:  Comparison  of  methods  on  a  128  x  128  cell  computational  grid  for 
the  vortex  flow  at  t  =  4. 
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(a)  Initial  circle 


(b)  64  x  64  grid  cells 


"Sir 


(c)  128  x  128  grid  cells 


(d)  256  x  256  grid  cells 


Figure  5:  HJ-WENO  particle  level  set  solutions  after  one  period  (t  =  8)  of 
vortex  flow. 
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(a)  Initial  circle 


(b)  64  x  64  grid  cells 


(c)  128  x  128  grid  cells 


(d)  256  x  256  grid  cells 


Figure  6:  Semi-Lagrangian  particle  level  set  solutions  after  one  period  (t  = 
8)  of  vortex  flow. 
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Figure  7:  A  uniform  grid  HJ-WENO  based  particle  level  set  three  dimen¬ 
sional  deformation  test.  Reprinted  from  [5]. 
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Figure  8:  An  octree  based  semi-Lagrangian  particle  level  set  three  dimen¬ 
sional  deformation  test. 
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4  Conclusions 


We  proposed  a  fast  semi-Lagrangian  based  particle  level  set  method  for 
the  accurate  capturing  of  passively  advected  interfaces.  By  utilizing  mass¬ 
less  marker  particles  nearby  the  interface  as  a  characteristic-based  diffusion 
control  mechanism,  we  were  able  to  forgo  the  use  of  higher  order  accu¬ 
rate  advection  schemes.  These  schemes  significantly  increase  the  computa¬ 
tion  time,  due  to  the  small  time  step  required  for  stability  of  the  scheme. 
Semi-Lagrangian  methods  do  not  suffer  from  this  constraint,  and  the  exces¬ 
sive  amount  of  numerical  diffusion  present  in  a  first  order  accurate  semi- 
Lagrangian  method  is  successfully  counteracted  by  the  particles.  Also,  se¬ 
lective  adaptive  mesh  refinement  for  small  scale  feature  resolution  is  compu¬ 
tationally  feasible  with  the  semi-Lagrangian  based  particle  level  set  method, 
unlike  other  more  complicated  Eulerian  based  advection  methods.  Moreover, 
the  computational  overhead  incurred  by  the  use  of  an  octree  grid  represen¬ 
tation  is  minimized  by  the  large  semi-Lagrangian  based  time  steps  allowed. 
Additional  optimization  is  gained  by  confining  the  level  set  update  to  a  nar¬ 
row  band  about  the  interface  and  the  use  of  a  fast  marching  method  to 
reinitialize  </>. 
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